/*
 *  sort.c
 *  Sandbox
 *
 *  Created by Joshua Lewis on 10/27/10.
 *   A parallel sorting algorithm (sample sort) using OpenMP.
 *
 */

#include "sort.h"

void psort(void *base, size_t num, size_t size, int (* comparator) (const void *, const void *))
{
  int threadNum = omp_get_max_threads(); // Number of threads OpenMP will spawn
	int sampleNum = 64; // Number of random samples to draw per thread
  int sortSize = num * 3; // Buffer size for each thread, N * 2 would occasionally cause problems
  int i, j, k;
  
  // Malloc chars, which we know are one byte, to store arbitrary data
	int *count = (int *)malloc(threadNum * sizeof(int));
	char *splitSample = (char *)malloc(threadNum * sampleNum * size);
	char *split = (char *)malloc((threadNum - 1) * size);
  // If the sampling does a bad job this will not be enough shared mem and weird things will happen
  char *threadSort = (char *)malloc(sortSize * size);
  
	for (j = 0; j < sampleNum * threadNum ; j++) // Randomly sample from base
    for(k = 0; k < size; k++)
      splitSample[j * size + k] = ((char *)base)[(rand() % num) * size + k];
	
	qsort(splitSample, sampleNum * threadNum, size, comparator); // Sort the samples
	
	for (j = 0; j < (threadNum - 1); j++) // Pick threadNum - 1 split points
    for(k = 0; k < size; k++)
      split[j * size + k] = splitSample[(j * sampleNum + sampleNum - 1) * size + k];
	
#pragma omp parallel private(i, j, k)
	{		
		int th_id = omp_get_thread_num();
    int th_index = th_id * (sortSize / threadNum) * size;

    // Search through the data and find the samples this thread needs to sort
    j = 0;
		for (i = 0; i < num; i++)
				if ((th_id == 0 || comparator(base + i * size, (void *)split + (th_id - 1) * size) > 0) &&
            (th_id == threadNum - 1 || comparator(base + i * size, (void *)split + th_id * size) < 1))
        {

          for(k = 0; k < size; k++)
            threadSort[th_index + j * size + k] = ((char *)base)[i * size + k];
          j++;
        }
		
		count[th_id] = j; // Number of samples found
                
		qsort(&threadSort[th_index], j, size, comparator); // Sort the samples
		
#pragma omp barrier
		
    j = 0;
		for (i = 0; i < th_id; i++)
			j += count[i];
  
    // Write the sorted samples back out to shared memory
		for (i = 0; i < count[th_id]; i++)
      for(k = 0; k < size; k++)
        ((char *)base)[(j + i) * size + k] = threadSort[th_index + i * size + k];
  }
	
	free(split);
	free(splitSample);
  free(count);
  free(threadSort);
  
	return;
}